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Abstract. In the context of Monte Carlo simulations, the analysis of the 
probability distribution P^(m) of the order parameter m, as obtained in 
simulation boxes of finite linear extension L, allows for an easy estimation of 
the location of the critical point and the critical exponents. For Ising-like systems 
without quenched disorder, P^(m) becomes scale invariant at the critical point, 
where it assumes a characteristic bimodal shape featuring two overlapping peaks. 
In particular, the ratio between the value of P^(m) at the peaks (Pi, max ) and the 
value at the minimum in-between (Pi jm i n ) becomes L-independent at criticality. 
However, for Ising-like systems with quenched random fields, we argue that instead 
AFl := In (Pl, mux / Pl, min) ^ L should be observed, where 9 > is the 
"violation of hyperscaling" exponent. Since 8 is substantially non-zero, the scaling 
of AFj, with system size should be easily detectable in simulations. For two fluid 
models with quenched disorder, AF^ versus L was measured, and the expected 
scaling was confirmed. This provides further evidence that fluids with quenched 
disorder belong to the universality class of the random-field Ising model. 



Fluids with quenched disorder 



2 



1. Introduction 

It was postulated by de Gennes that liquid-gas type transitions in fluids confined 
to quenched random porous media belong to the universality class of the random 
field Ising model (RFIM) [l]. The unambi guous experimental verification of this 
conjecture remains elusive to this day [2], but there is growing numerical evidence 
from computer simulations [3j|5] of quenched-annealed mixtures |6]. In quenched- 
annealed mixtures, a fluid of mobile particles is confined to a configuration of immobile 
(quenched) particles. The quenched particles induce a random field provided (1) their 
positions are sufficiently random and (2) they display a preferred affinity to one of 
the phases formed by the mobile particles. When either one of these conditions is not 
met, the conjecture of de Gennes does not apply |2j[7]. 

Computer simulations of quenched-annealed mixtures remain complicated. One 
problem is that for systems exhibiting quenched disorder an additional average over 
many samples drawn from the distribution characterizing the disorder needs to be 
taken. In addition to the thermal averaging, a disorder averaging over the different 
samples must be taken, and hence the computational effort is orders of magnitudes 
larger. The convergence of disorder averaged quantities with the number of samples 
is typically slow (random field systems at criticality do not self average [8}fl3]) and so 
the disorder averaging must comprise many samples. A second problem is that the 
critical exponents of the three-dimensional RFIM are not known very precisely. In 



particular, estimates for the correlation length exponent v range from 1.1 to 2.25 14 
Due to the large uncertainty in v, standard finite size scaling (FSS) of simulation 
data is problematic. The only critical exponent of the RFIM on which there is some 
consensus, is the "violation of hyperscaling" exponent 9. In pure systems, meaning 
systems without quenched disorder, it holds that 6 = 0, but in the RFIM 9 ~ 1.5, 
i.e. distinctly non-zero. In order to detect RFIM universality, it follows that FSS of a 
quantity "somehow involving" the violation of hyperscaling exponent 9 is a promising 
approach. In this paper, we present one such strategy, based on the free energy cost 
of interface formation. In the RFIM at its critical point, the barrier diverges with the 
system size as 

AF L cxi e , (1) 
where L is the lateral extension of the simulation box |5|. In the pure Ising model 



b ox [5], 



where 9 = 0, AFl is L-independent at criticality |5||15| 16 . Since 9 is large in the 
RFIM, the divergence of AFl should be easily detectable in random field systems. 

The outline of this paper is as follows. We first describe how the scaling of the free 
energy barrier at criticality can be derived and exploited in the context of FSS. We then 
illustrate the approach for the RFIM, the Widom-Rowlinson (WR) model [17] with 
quenched obstacles, and the Asakura-Oosawa (AO) model [l8|[T9] of colloid-polymer 
mixtures inside a porous medium. In all cases 9 > is observed, but the precise value 
varies. For the RFIM and WR mixture 9 ~ 1.5 is obtained, while colloid-polymer 
mixtures yield a somewhat lower value 9 ~ 1.0; possible origins for this discrepancy 
are discussed in Section 2] 
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2. Finite size scaling of the interface free energy 

2.1. theoretical background 

We assume that the RFIM in d = 3 dimensions undergoes a continuous phase 
transition, at critical temperature T c , from a disordered phase at high temperature, 
to an ordered phase with finite magnetization at low temperature. The existence of a 
transition at nonzero temperature was controversial until a proof for the existence of a 



spontaneous magnetization settled this issue 20 ; rigorous results on the order of the 
transition remain elusive, but numerical studies [21 22 indicate that the transition 
is continuous. Following convention, we define the relative distance from the critical 
point as 

t := T/T c - 1. (2) 

In the vicinity of the critical point, t = 0, the correlation length diverges as a power 
law 

£oci*r", (3) 

while in the ordered phase a finite interfacial tension (excess free energy per unit of 
interface) develops 

a cx \t\ 2 - a -» (t < 0), (4) 

where a is the critical exponent of the specific heat. Near the critical point, the 
correlation length is the only relevant length scale, and so a simple dimensional 
argument implies that the total interface free energy must scale as 

AF$ CX CT^" 1 CX |f|2-a-r/^-l K ^ a +du-2)/u^ (g) 

where in the last step t was eliminated using Eq. (§. If we now use the FSS "Ansatz" 
£ cx L, the above equation reduces to 

AF L (xlK^)/^ ( 6 ) 

which describes the scaling of the interface excess free energy in a finite simulation 
box of lateral extension L in the vicinity of T c . This equation applies to both the 
pure Ising model and the RFIM. The hallmark of RFIM universality is the modified 
hyperscaling relation [23T|25] 



2-a = v(d-0), (7) 

which upon substitution in Eq. ([6| yields 

AFl cx L e , (8) 

and so we have derived Eq. ([lj. In the pure Ising model = 0, in which case the 
barrier does not depend on L at the critical point (as is well known |15| ). In the 
RFIM ~ 1.5, implying a strong divergence with system size. 
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2.2. order parameter distributions 

We now explain how the barrier AFl is measured in computer simulations of liquid- 
gas type transitions. For simplicity, we first consider a single component fluid in the 
absence of quenched disorder. The key quantity is the order parameter distribution 
(OPD) 

Pl(p\T,/j.) (single component fluid), (9) 

defined as the probability to observe a state with particle density p. The OPD is 
measured in the grand canonical ensemble, at fixed temperature T and chemical 



potential p, using a periodic L x L x L simulation box 26 27 . Below the critical 
temperature T c and at the coexistence chemical potential p cx , the OPD becomes 
bimodal; the peak at low (high) density corresponds to the gas (liquid) phase. The 
coexistence chemical potential is obtained by maximizing the derivative of the average 
density with respect to p 

Mcx : d(p)/dp -> max, (10) 

with (p) = J p Pl(p\T, p) dp. Note that other choices to define p cx are possible 



also [28] , most notably the "equal-weight" rule [29], but Eq. (10 1 has the advantage 
that it remains well-defined also when the peaks in the OPD overlap. The free energy 
barrier AF^ is encoded in the logarithm of the OPD: it corresponds to the average 
peak height of \uPl(p), measured from the minimum "in-between" the peaks [30] , 

In athermal binary mixtures with mobile species A and B (but still without 
disorder), the analogue of Eq. ([£]) becomes 

Pl(pa\za,zb) (athermal binary mixture), (11) 

where pa is the density of the A particles, and where za and zb arc the fugacities of the 
A and B particles, respectively. By athermal we mean that the particle interactions are 
either hard-core or ideal, as is the case in the WR and AO models. Those mixtures 
correspond to a single component fluid if one identifies p f-> pA, T <-> 1/zg, and 
p «-» lnz^. The transition between a "gas phase" with low pA and a "liquid phase" 
with high pa is now driven by the fugacities. It occurs on the coexistence curve 
za = Zcx(zb), provided the fugacity of the B particles exceeds the critical fugacity 



zb > z c . On the coexistence curve, which again can be found using Eq. ( 10 1, the OPD 
becomes bimodal and allows the extraction of a free energy barrier AF^, see Fig.JTJa). 

We now consider an athermal binary mixture in the presence of quenched disorder. 
We thus have mobile particles A and B as before, but also quenched obstacles which 
are distributed at random locations at the start of each simulation after which they 
remain frozen. The A and B particles then diffuse through the quenched environment. 
From the set of all possible configurations of quenched disorder, we consider a finite 
number of N configurations, and for each of them the OPD is measured. Instead of a 



single OPD as in Eq. (Ill, we thus have a set of distributions 



Pla{pa\zb,z a ), (12) 

where i = 1, . . . , N labels the different configurations. Even though we are interested 
in the OPDs for different fugacities, it is sufficient to simulate only a few combinations: 
having simulated a system at some (za,zb), the OPDs for nearby fugacities can be 
obtained via histogram reweighting 1 3 1 . To facilitate an efficient storage of the data, 
the approach of |3l^j is used. Meaningful results around the critical point require 
N ~ 10 3-4 , so the computational effort is much larger compared to the pure systems. 
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Figure 1. (a) Example OPD for the WR mixture (with quen che d obstacles) 
obtained for zg = 1.2, L = 12, and z& chosen according to Eq. Note that 

the natural logarithm of the distribution is shown. The barrier AFl i corresponds 
to the average peak height (vertical arrow), (b) Example of two "problematic" 
OPDs, again for zg = 1.2, L = 12, and Za tuned according to Eq. jlQ) . These 
distributions feature a peak at intermediate density, which does not correspond 
to a gas or liquid phase. The intermediate peak either replaces the gas or liquid 
peak, in which case the resulting distribution remains bimodal (solid curve), or it 
appears as a third peak (dashed curve). A free energy barrier corresponding to 
coexisting gas and liquid phases cannot be extracted from these distributions. 



2.3. extracting the free energy barrier 

To obtain the quenched-averaged free energy barrier, we extract the free energy barrier 
AFr^ for each configuration separately, as in Fig. [Tja), and then average these values. 



To be precise: for a given zb, each OPD is tuned to coexistence via Eq. (10) 
which implies that z cx varies between disorder configurations [3j[5]. By tuning the 
distributions separately, we ensure that the majority of them become bimodal, such 
that a barrier can indeed be "read-off". From the individual barriers, the quenched- 
averaged barrier is obtained MS i\ simple average: 

N 

AF l = (1/N)J2&Fl, 1 - (13) 

i=l 

The finite number of disorder configurations used in calculating the barrier gives rise 
to a statistical deviation 

1 N 

which represents the expected deviation of our result from the result we were to obtain 
for N — > oo. 

In practice, obtaining the barrier AFr^ can be problematic. For some disorder 
configurations, the OPD does not feature a clear liquid and gas peak under the criterion 



of Eq. ( 10 1. In these cases, an intermediate peak is seen, at a density roughly between 



that of the liquid and gas. This extra peak can occur as a replacement for cither 
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Figure 2. Illustration of how the "problematic" OPDs in Fig. [TJb) might 
arise. We provide data for the WR model (particle species A and B) with a 
configuration of quenched obstacles explicitly tailored to feature a peak in the 
OPD at intermediate density. Instead of distributing the obstacles randomly, 
they are placed on a regular grid; the obstacles in the lower half of the grid favor 
A particles, those in the upper half B particles. This configuration yields three 
stable states, depicted in the snapshots of the upper row, which were obtained for 
zg = 1.3 and L = 12. For clarity, only the A particles and obstacles are shown. 
The first state is a homogeneous gas phase, characterized by a low density of A 
particles (left frame). The second is an inhomogeneous state, with a high (low) 
density of A particles in the lower (upper) half of the box (center frame). The 
third is a homogeneous liquid phase, characterized by a high density of A particles 
(right frame). The graph shows the corresponding OPD at three values of za- 
For za = zg, a single peak is observed at intermediate density (solid curve); for 
za < zg, a bimodal distribution is obtained with the second peak appearing at 
low density (dashed curve); for Za > zg, a bimodal distribution is also obtained 
but with the second peak appearing at high density (dotted curve). 



the liquid or gas peak, or as an additional third peak (Fig. [TJb)) . The existence 
of these "problematic" distributions is due to cavities formed by the obstacles. In 
some disorder configurations, the obstacles create large regions with a local preference 
for either A or B particles. One can easily show that such regions induce an 
intermediate peak (Fig. |2|. We emphasize that around the critical point such 
problematic distributions are scarce in the full ensemble of disorder configurations 
(although modified hyperscaling does allow a finite fraction of them pH). 
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2.4- analysis of finite- size effects: scaling plots 

Near the critical point, we expect AFl to scale according to Eq. ([lj. Following 
standard FSS practice (32 33 , we test Eq. ([T| using scaling plots. That is, we plot 
Ul '■= L~ 9 AFl versus r := tL 1 ^, t given by Eq. (I2J. Then, if the correct values for 
9, is, and T c are used, the curves j/l(t) for different system sizes L should collapse on 
top of each other. The practical difficulty with this method is assessing the quality 
of the collapse: seemingly good data collapses are obtained over a range of parameter 
values, and simple "eye gauging" becomes unreliable. To obtain the parameters of best 
collapse, as well as a measure of the reliability of the results, a numerical procedure 
was therefore used. 

Given a candidate tuple (9, v, T c ), Eq. (14 1 directly yields the statistical 
uncertainty in Ul{ t ), namely 8ul(t) = L~ 6 Ul- We average 6ijl(t) over the different 
system sizes to obtain the typical statistical uncertainty u(t) of the curves at value r. 
For a given t, the curves i/l(t) for the different system sizes spread around some 
average value with a root-mean-square width ct(t). We then compare <t(t) to the 
typical statistical uncertainty u(t) via the ratio cr(r) /u{t). This ratio measures 
how strong the curves differ compared to the random scatter induced by the finite 
number N of obstacle configurations. To assess the quality of the collapse, we then 
compute the average of the ratio over some region of r around the critical point 

R(0, u, T c ) := r 44 dr. (15) 

Tl - T J TQ u{t) 

Values R < 1 are interpreted as fully consistent with a perfect collapse of the curves, 
while R > 1 is considered less-consistent. The integration range is chosen as large as 
possible, but small enough to stay within the critical region. We used r = — 0.2- f CP 1 /" 
and n = 0.2 • fCr 1 /", which assumes a typical system size L ~ 10, and restricts the 
temperature range to \t\ < 0.2. 



3. Results 

3.1. RFIM 

To provide a benchmark point for our analysis, we first apply our scaling test to 
data for the RFIM model. These data were taken from our previous work (5], 
and obtained for the RFIM in three dimensions using a random field drawn from 
a Gaussian distribution. The total number of disorder configurations per system size 
equals N ~ 10 4 . For system sizes L = 8, 10, 14 and 16, we examine R(9, v, T c ) 
in the plane of critical exponents (9, v). For each point in the plane, the critical 
temperature T c is tuned such that R{9, v,T c ) is minimized; the value at the minimum 
is denoted R m in(9, v). The results are collected in Fig.[3]as contours. The region where 
-Rmin(0, v) < 1, i.e. where the collapse is fully consistent, agrees with the expected value 
9 ~ 1.5. Also of interest is that the region of good collapse features a tail extending 
toward the pure 3D Ising values (9 — 0, v ~ 0.63) [34]. The cause of this tail is not 
entirely clear, but we suspect it is due to crossover effects (35] . However, at the pure 
Ising values, the deviation between the curves j/z, exceeds the statistical fluctuation 
due to the finite number of disorder configurations by a factor of more than three. 
We conclude that the pure Ising exponents are effectively excluded by our analysis. 
Note that 9 = d — 1 = 2, i.e. the "exponent" one would observe below the critical 
temperature (where the transition is first-order) is also excluded. 
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Figure 3. Scaling analysis for the RFIM. Shown are contours of R m i n (0, u) in 
the plane of critical exponents (0,v); dots mark the critical temperature T c at 
a few selected points, the best estimate being T c ~ 3.4. The plot shows a clear 
preference for the expected RFIM value 9 ~ 1.5, while 9 = of the pure Ising 
model, as well as 9 = 2 of a first-order transition, are effectively excluded. 




Figure 4. The analogue of Fig. [3] for the WR model with quenched obstacles; 
dots indicate the critical fugacity z c at some selected points in the (9, v) plane. 
The global minimum of R m i n (9, u) is at z c ~ 1.42, which thus serves as our best- 
estimate for the critical fugacity. As for the RFIM model, a clear preference for 
9 ~ 1.5 is found, as well as the exclusion of 9 = and 9 = 2. 
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3.2. WR model with quenched obstacles 

The WR model [17] consists of unit-diameter spheres, of species A and B. The only 
interaction is a hard-core repulsion between particles of the opposite species. We 
consider the WR model in the presence of quenched obstacles; the spatial dimension 
is d = 3. The disorder configurations are generated by placing 0.02 x L d particles 
of species A, and the same number of species B, into a periodic cube of edge L 
(non-integer numbers are appropriately rounded to the next integer at random). The 
obstacles are inserted at random positions, irrespective of overlap between them, after 
which they remain quenched. Next, mobile particles are introduced and these interact 
with the quenched particles following the usual WR rule: quenched A-obstacles (B- 
obstacles) have a hard-core interaction with mobile B (A) particles. In this work, 
system sizes L = 7, 8, . . . , 12 are used. For each system size, OPDs for N ~ 10 4 
disorder configurations were generated at fugacities za = zb = 1.4, which is the 
approximate location of the critical point as anticipated from our previous simulations 
[5]; all details regarding the simulation method can also be found in that reference. 

To obtain AFl via Eq. (T3| over a range of zb requires extracting the free energy 
barrier from several million(!) distributions, which obviously cannot be done by hand. 
To deal with the "problematic" distributions of Fig. [TJb) automatically, a numerical 
filter was therefore applied: if one of the two dominating peaks in the OPD is centered 
in the intermediate density range 0.20zb < PA < 0.75zb, the distribution is considered 
to not feature a distinct liquid and gas peak. For these distributions we set AFlj = 0. 
Since for the values of zb considered here no more than ~ 0.1% of the distributions 
failed this criterion, our results should not significantly depend on the filter. The 
results for R m i n (9, v), as well as some estimates for the critical fugacity z c obtained 
from minimizing R, are shown in Fig. [4] Note that, for the WR model, the analogue 
of Eq. |2| becomes t := z c /zb — 1. As in the RFIM, Fig. [IJfeatures a clear preference 
for 9 ~ 1.5, and the competing values 9 = and 6 = 2 of, respectively, the pure 
Ising model and a first-order transition, are effectively excluded. We interpret this 
finding as a strong sign that the WR model with quenched obstacles indeed belongs 
to the universality class of the RFIM. Note, however, that we cannot provide high- 
precision estimates of the critical exponents, nor of the critical "inverse temperature" 
z c . There is a substantial range of possible values (9, v, z c ) in the region R m i n (9, v) < 2, 
i.e. the region where our results for the WR model and RFIM become quantitatively 
compatible. 

3.3. Colloid-polymer mixtures in porous media 

We now consider a colloid-polymer mixture confined to a porous medium. The 
colloids and polymers are unit-diameter spheres, with only a hard-core repulsion 
between colloid-colloid and colloid-polymer pairs (this is just the AO model fl8)[l9"] 
with equally sized colloids and polymers). As porous medium we use a quenched 
configuration of polymers, which are distributed randomly in a periodic cube of edge 
L at the start of each simulation run (the average packing fraction of the quenched 
configuration r]Q = 0.05). The OPD for the AO model is the probability distribution 
PL,i(Vco\\ z co\j lp) |27l, where i denotes the quenched polymer configuration for which 
the OPD was measured, rj co \ the colloid packing fraction, z co \ the colloid fugacity, 
and ?7p the polymer reservoir packing fraction. Note that, for the AO model, rf p is 
just the polymer fugacity multiplied by the volume of a single polymer, and thus 
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Figure 5. Best-estimates of (t9, u, r?p cr ) versus the scaling variable r for the AO 
model inside a porous medium; see details in text. 



plays the role of inverse temperature. Consequently, the analogue of Eq. ([2| becomes 
t := ijp/r/p C1 . — 1. One of us (RV) recently studied the above model using the same 
type of quenched disorder. We now extend the analysis of that previous work [3j to 
the scaling of the free energy barrier. 

One practical problem is that the number of disorder configurations in our AO 
data set is "only" N ~ 10 3 , i.e. one order of magnitude below the value used in the 
analysis of the RFIM and WR model. Also the quality of individual distributions 
was worse compared to that of the latter two models. A further complication is that 
the AO model is asymmetric: the coexistence colloid fugacity is not trivially related 
to 77p. In contrast, for the WR model, it holds that za = Zb at the critical point 
due to symmetry. Because of all these problems, the construction of contour plots, 
such as Fig. [3] and Fig. |4j was not feasible for the AO model, and so a modified 
analysis is presented. To this end, we extract the barrier AFt, from the quenched- 
averaged OPD, following the method outlined in the Appendix of Ref. [5]. This 
procedure yields a result similar to Eq. (13), but it is less susceptible to statistical 
uncertainties in individual distributions. Since the barrier is extracted from the 
quenched-averaged OPD, we loose information about the fluctuations between disorder 



configurations. In particular: Eq. (14 1 can no longer be used to calculate the typical 
statistical uncertainty u(t) defined in Section 2.4 Hence, in the construction of the 



scaling plot yL( T ), we cannot compare the deviation between the curves cr(r) to the 
statistical uncertainty u(t). Only best-estimates can be provided, which are obtained 
by minimizing the relative deviations between the curves. That is, for a given r 
we minimize the value a(r)/y{r), where y(r) is the average value of the j/l(t). For 
—0.05 < r < 0.05, and system sizes L = 10, 11, 12 we computed these best-estimates: 
their variation with r is shown in Fig. [5] If the collapse was perfect, the best-estimates 
should not depend on r. Regarding rf p cr and 6 this is confirmed, and from the average 
of the curves we conclude that rf v cr ~ 1.3 and 9 ~ 1.0 for the AO model. However, it 
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is clear from Fig. [5] that no meaningful estimate of v can be provided. We emphasize 
that r/p cr ~ 1.3 deviates by about 10% from previous estimates. The reason is that, 
in Ref. pH, we assumed i^rfim = 1-1. However, from the analysis of the RFIM (Fig. [3]) 
it transpires that vrfim may well be different. This is furthermore corroborated by 
the literature |14|, where large variations in z^rfim are also reported. 



4. Conclusions 

The aim of this paper was to test the relation AFl oc L in fluids with quenched 
disorder. This relation applies when the universality class is that of the RFIM j5]. 
For the WR and the AO model with quenched obstacles, our data confirm that AFt, 
diverges as a power law with system size at the critical point. This is an important 
qualitative indication that the universality class is no longer that of the pure Ising 
model, since then AFt, would be constant at criticality. For the WR model, 9 is 
consistent with the RFIM value #rfim ~ 1-5. For the AO model, 9 ~ 1.0 is found, 
which is somewhat below the RFIM value. Possible reasons for this discrepancy are 



(1) crossover effects 35 (in this case from pure Ising to RFIM universality), (2) 
corrections to scaling due to asymmetry, and (3) insufficient disorder averaging. We 
believe that the crossover scenario is the most likely explanation, since 9 ~ 1.0 is "in- 
between" the pure Ising value 6 — and the RFIM value. Of course, all these problems 
could be solved by simulating larger systems and more disorder configurations, but 
the cost in CPU time is still prohibitively large. Our data also make clear that high- 
precision estimates of critical exponents and temperatures in random-field systems 
remain difficult to obtain in simulations. The data for the RFIM and the WR model 
reveal large variations, see Fig. [3] and Fig. |4j in particular for the correlation length 
critical exponent v. 
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